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A METHOD FOR SOLVING HYPERVELOCITY IMPACT PROBLEMS 
WITH CONSIDERATION OF FREE-SURFACE DISCONTINUITIES 


By Richard E. Turner 
Langley Research Center 

SUMMARY 

A new method of tracking free surfaces has been developed and used to apply 
boundary conditions along free-surface discontinuities occurring in the motions of com- 
pressible media. The stability of the finite-difference equations used to track those free 
surfaces has been investigated under simplified conditions and the method was found to 
be stable. 

Also, finite -difference approximations to the continuity equation and general trans- 
port equation (which governs the motions of continuous media, in cylindrical coordinates) 
have been developed in which convective terms are simplified and motion constants (mass, 
momentum, and internal energy) are conserved in transit. 

This system of equations has been applied to the solution of a hypervelocity- impact 
problem in which an aluminum cylinder with a 10-centimeter base diameter, 8-centimeter 
top diameter, and 6-centimeter height impacts a 2-centimeter-thick target with a relative 
velocity of 20 kilometers per second. This solution is compared with an independent 
solution. 


INTRODUCTION 

The acquisition of the scientific knowledge necessary to insure that future space- 
craft have adequate (but not excessive) meteoroid shielding will require at least three 
distinct research approaches. Space-flight experiments are necessary to obtain a lim- 
ited knowledge of the damage potential of the meteoroid environment. Laboratory impact 
experiments are necessary to infer the meteoroid environment from the flight meteoroid- 
damage experiments. Analytical studies of the meteoroid impact-damage phenomena can 
then be used to complement a limited capability for laboratory impact simulation so that 
a better knowledge of the damage potential of the meteoroid environment can be obtained 
and protection against it can be provided. 

A variety of digital computer programs have been developed expressly for the pur- 
pose of studying these hypervelocity impact problems. Some programs treat the working 



media as a compressible fluid, some as elastoplastic solids, and some as perfectly 
plastic solids. The aim of these calculations is to find crater sizes, deformations on 
the rear face of finite -thickness targets, and deformations of secondary targets in the 
case of multiwalled structures. The coordinate viewpoint may be Eulerian, Lagrangian, 
or Eulerian-Lagrangian mixtures. (For some examples, see ref. 1.) 

In programs with the Eulerian viewpoint, the free surfaces of the projectile-target 
system are treated by diffusion. Discontinuities at the free surfaces are smeared out 
over several nodal distances. Reference 2 is an example of such a treatment. The 
"particle in cell" approach (a dual use of Eulerian and Lagrangian viewpoints) also treats 
the free surfaces by diffusion. Some cases calculated by this approach are presented in 
reference 3. The method of marker particles is sometimes used to locate free surfaces; 
however, the marker particles are only visual aids and do not affect the free-surface 
treatment. 

The treatment of free surfaces by diffusion has some disadvantages as well as 
advantages. Two strong advantages are simplicity and the ability to approximate com- 
plex free-surface geometries. Two strong disadvantages are that boundary conditions 
cannot be applied accurately and the position of the free surfaces cannot be defined 
accurately. 

A method for tracking free surfaces is given here which can be used to apply bound- 
ary conditions along the free surface of a system provided the slope of the free surface 
is finite everywhere. A mass transport equation and a general transport equation are 
derived in the Eulerian frame of reference from an intuitive basis. These transport 
equations represent a simple treatment of the convectivelike terms in a continuum and 
yield smooth space solutions to the finite -difference equations. 

The model used to represent the target -projectile medium is a compressible non- 
conducting fluid with a spherical stress tensor. (See ref. 4.) 

SYMBOLS 

A,B time -dependent coefficients in a Fourier series, centimeters/second 

a,b sgn(-u), sgn(-w), dimensionless 

C time -dependent coefficient in a Fourier series, centimeters 

E internal energy per unit mass, ergs/gram 

f,g equations of state 
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H 


damping constant, dimensionless 


M 

P 

Q 

q 

r 

S(p) 

S(pE) 

S(pQ) 

S(pu),S(pw) 

s 

t 

U 

u 

V 

W 

w 

X 

z 


amplification matrix appearing in numerical stability studies, dimensionless 
pure hydrostatic pressure, dynes/ centimeter 2 


arbitrary quantity being transported, dimensions arbitrary 
artificial viscosity used to generate entropy, dynes/centimeter2 


radius, centimeters 

density source, grams/centimeterS 


energy source. 


ergs/centimeter ^ 


second 


source for arbitrary quantity Q being transported, 
(grams/ centimeter 3) (Q dimensions) 
second 


sources for momentum, g ram s/centimeter2 

second2 

displacement along a free surface, centimeters 
time, seconds 

radial-velocity disturbance, centimeters/second 
radial velocity, centimeters/second 
velocity, centimeter s/second 

longitudinal-velocity disturbance, centimeters/second 

«► 

longitudinal velocity, centimeters/second 
specified radial position, centimeters 
longitudinal coordinate, centimeters 
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|3 space translation, centimeters 

y Fourier exponent used in numerical stability study, centimeter - * 

6 elemental volume in cylindrical coordinates, centimeters^ 

0 angular variable in cylindrical coordinates, radians 

6 angle between the free surface and a constant-z parametric line, radians 

X eigenvalues of a square matrix, dimensionless 

£ specified radial position, centimeters 

p density, grams/centimeter^ 

a pure tensile stress, dynes/centimeter^ 

t time increment, seconds 

<f> free-surface-coordinate disturbance, centimeters 

0 free -surface coordinate, centimeters 

0 state variable for target-projectile system, arbitrary 

Subscripts: 

B evaluated along a free surface 

j radial coordinate jAr 

k longitudinal coordinate kAz 

Superscript: 

1 time It 


Notations: 



square matrix 
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column vector 

Euclidean norm, dimensions arbitrary 


The operator d( )/dt denotes time derivative following the particle (comoving), 
A bar over a symbol denotes a vector. 


ANALYSIS 


The problem to be solved concerns a solid circular cylinder impacting a solid tar- 
get of finite thickness. The pressure generated by the impact is assumed to be suffi- 
ciently high to allow yield strengths to be neglected. The resulting medium can be 
approximated by a compressible fluid. 


Differential Equations of Motion For Compressible Nonconducting Fluids 

With Pure Hydrostatic Stress 

The differential equations of motion for compressible fluids with pure hydrostatic 
stress are (from ref. 4) the momentum equation 


the energy equation 


and the continuity equation 



where p is mass density, V is velocity, a is stress, E is internal energy per unit 
mass, t is time, d( )/dt is the comoving time derivative of ( ), and V is the classi- 
cal vector del operator. The operator d( )/dt is given by 

= 2D. + v . v( ) (4) 

dt °t 

where 9( )/8t is the time derivative of ( ) for fixed-space coordinates. The algebraic 
equations needed to complete the description of compressible fluids with pure hydrostatic 
stress are the constitutive equation • * 

a = -(p + q) (5) 
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where p and q are nonnegative, and the equation of state 


P = f(E ,p) (6) 

In these equations p is hydrostatic pressure and q is an artificial viscosity term 
which is introduced for entropy production in such things as shock waves. (See ref. 4.) 


Interior Equations of Motion in Finite -Difference Form 

If r is the radial coordinate, z is the longitudinal coordinate, u is the radial 
velocity component, w is the longitudinal velocity component, A denotes a coordinate 
increment, the superscript l refers to time It, the subscript j refers to the radial 
coordinate jAr, k refers to the longitudinal coordinate kAz, and the fraction 1/2 in 
subscripts refers to values averaged from interval end points, then from appendix A the 
equations of motion for the medium can be written as the transport equation 


J 


j>k at 


= -a 


p j+(a/2),k u j+(a/2),k 


Q 


j. + a ? k 






Ar 


-bpi /os W Z 


Q' 


.i'k+b 


j,k+(b/2) j ,k+(b/2) A z 


Q* 

_ILiS 


Q' 


l 




p j+(a/2),k u j+(a/2),k I+a ’j J ’ k + s <pQ)^ k 


- r 




at 


at 


( 7 ) 


where Q assumes the role of E, u, and w, and where a is sgn(-u) and b is 
sgn(-w). The source equations are 


s <H,k 


Ar 


( 8 ) 


S(p W )i , = j, Ml/2) - g i,k-(l/2) 
j Az 


(9) 






dt 


( 10 ) 
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along with the continuity equation 


!ik + au* , /9 x^ +a > k ~ P ^ k + bw* f. /9 /l k+b ~ P l k + u^ , /9 v . P l^ ? k " P U 

at ]+(a/2),k Ar ],k+(b/2) A z ]+(a/2),k 2 rj 


_ i u U/2),fc ~ u )-(l/2),k l jb(l/2)-^-(l/2) i u W/2),k tu j-(l/2),k 

P 3 ,k Ar P ],k Az p ],k 2r . 

( 11 ) 


and 


dp! ir ^p. 7 

— “= = — 1*== + a u. / /os , 
dt 8t j+(a/2),k 


W p i,k + b1 j 


l l 

p j,k+b " P j,k 


Ar 


j,k+(b/2) az 


+ n l , /9 s . P 3 +a /k-L P lk 

j+(a/2) ? k 2rj 


( 12 ) 


From reference 5, the equation of state for aluminum, valid to about 
5 x 1C)12 dynes/centimeter 2 , is 



if p ^ 0. If the equation yields negative values, p is taken as zero. 

The damping term q, which is introduced to render the numerical solution stable 
and give structure to the shock wave (see ref. 4) is 

H(Ar 2 + Az 2 ) | ^ | dp 

q " P "dt” dF (ld) 

where H is a constant which is chosen by trial and error. The value for H used here 

was 5 for ^ £ 0; otherwise H is zero, 
dt 

Boundary Considerations 

Equations (7) to (12) can be solved in a straightforward manner if the fluid medium 
occupies a simple cylindrical region. In the problem considered here the fluid medium 
occupies a region which has little resemblance to a cylinder. The approach often used to 
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circumvent this difficulty is to surround the interior region with a cylindrical region. 

The region between the problem boundary and cylindrical boundary is then assumed to be 
a vacuum. (See fig. 1.) As time increases, the projectile moves into the target region 
and both regions become distorted and spill into the empty regions. In taking differ- 
ences along the free surfaces, the discontinuities in the system properties are ignored. 
Here, along the free surfaces, the requirement for smooth solutions and the simulta- 
neous requirement for free-surface definition meet in conflict. Smooth solutions require 
that short-wavelength space variations be attenuated while good free-surface definition 
requires that short-wavelength space variations be maintained. This incompatibility 
between smoothness and free-surface definition can be eliminated by tracking the free sur- 
faces and taking the boundary discontinuities into considerations near the free surfaces. 


Free-Surface Equations of Motion 

Consider the surface of revolution shown in figure 2, which separates the interior 
from the void region. This surface can be tracked by tracking the intersection of the sur- 
face of revolution with the vertical parametric lines of constant radial coordinate. Let 
0 be the distance along a vertical parametric line to the surface of revolution, so that 

0 = <M r jt) (14) 

It is seen from figure 3 that 

0(x, t+r) = 0(|, t) + wr (15) 

Let A 0 be 0(x, t+r) - 0(x,t ); then equation (15) is equivalent to 

A0 = w ir (16) 

and since x - £ equals ur, equation (16) can be divided by r to get 


A0 _ 0(x,t) - 0(g, t) 

r X - 4 


(17) 


Now make the approximation that r tends to zero, and transform equation (17) to 
the index notation to get 


90 1 

“at 




1,B ' aU j,B 


,1 ,1 

^-i+a " 

Ar 


(18) 


Since 0^ will not, in general, coincide with any of the system's interior points, 
equation (18) introduces two quantities to keep track of: u^ R and v/'. Equations for 

these quantities can be developed from equations (1) to (6). Two approaches will be given 

for computing u^ „ and w^ D . 

lit* 
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Longitudinal, coordinate, 


Void 


Target 


Void 


\ 


Projectile 


Radial coordinate, r 


Figure L- Typical initial geometry of an impact problem. 



Radial coordinate, r 


Figure 2 - Typical surface of revolution intersecting the parametric lines of constant radius 

in cylindrical coordinates. 



Figure 3.~ Motion of surface of revolution relative to the vertical lines of constant radius. 


First approach .- Along the free surface, a is zero. From equation (5), along the 
free surface, 

a = -(p + q) = 0 

Thus, since p and q are nonnegative, 

p(Boundary, t) = q(Boundary, t) = 0 

Equation (2) shows that 

^(Boundary, t) = 0 
dt 

so that 


E(Boundary, t) = 0 


The equation of state, equation (6), can be inverted to give 

P = g(P,E) 


Since p = E = 0 along the free surfaces, 

p(Boundary, t) = g(0,0) = Constant. 

Equation (1) can now be used to find the two other equations needed for 

Consider figure 3 for a particle on the free surface, moving with 
The identity 


(19) 

u* R and w* R . 
the free surface. 
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V(x, t +r) - V(£,t) = [v(x, t+r) - V(x,t)] + [v(x,t) - V(|,tj] 


can be written, and the following equation can be derived in a manner similar to that in 
appendix B: 


'(¥' t+ 5 ) 


traiectom 


t Xi t+ ?) 


<J time V(x,t) - Vi 
r x - £ 


.x-1 


As in appendix B, the approximation is made that r tends to zero. Equation (21) is 
changed to index notation to get 

dv* av* V* -¥ 

!!k = AB tau ( D V fe B J ;B (22 

dt 8t ],B Ar V 

This expression is the comoving derivative of a point on the surface, and thus z is not 
an independent variable, but z = 0(r,t). 

When equation (22) is combined with equation (1), there follows 

dll\ R ; R “ ^ r. 1 do 1 . R 

— = -au^ tj 12 ® + _± -ii§ (23 

at J,B Ar n l 8r v 

P ],B 


Hb, 1 8<T j.B 

at J,B Ar + D l dz 

P i,B 

da \ B dQ \ B 

Values for ■■ and — can be found by noticing that 


ao* „ do 1 . 


_^ = __LB si ne + _-i£cose = o 
as 8z ar 


where s is a displacement along the free surface and 6 is the angle measured from 
the horizontal to the free surface in figure 4. The partial derivative of o\ „ with 

respect to z is 


9z kAz - 0^ 
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where 


Be 1 . 


LB 


0S 


a\ , is evaluated at the closest interior point, 
j 

and letting 


tan 0 = 


A 


i+l 

2Ar 


A-! 


Now solving for 


fig 

dr 


from 


yields 




Figure 4- Free-surface coordinates in relation to the fluid interior nodal points. 


12 



Second approach .- In the second approach equations (23) and (24) are still valid, 

Bp\ B _j 7 

but an equation for — is obtained from equation (22) by replacing V\ „ with p. „ 

to get 


°p\ 




dt 


at 


i,B 


Ar 


(27) 


d p\ 


An equation for — is obtained by noticing that the left-hand side of equation (A4) is 
dt 

dp i B 

equal to , and taking the limit of the right-hand side as Ar and Az tend to zero. 

The result is 


Now let 


dt *J> B l 8r az rj 


(28) 


^■B - U ll.B - U 1-1,B 
as As 


(29) 


9z <$k - kAz 

J 


(30) 


and 


9z 0? - kAz 


d>{ i - d& 1 

tan 6 = 3 +1 ■ -Slrl 

2Ar 


(31) 


(32) 


Then, since 


8u 


l 


8u- b R 

_li2 = — LE. cos 8 + — kE sin q 
as ar az 


(33) 


VU • T> 

solving for — from equation (33) and substituting from equations (29) to (32) gives 
ar 
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(34) 


"“j-B - “i+l.B ' (^,B- u j,k)(»jri - *i-l) 


2Ar 


kAzJ2Ar 


a P \ B 

Equations (27), (28), (30), and (34) define — £ — in terras of known quantities. 
' 8t 


Application of Boundary Conditions 

Figure 5 shows a target -projectile system during the time of crater formation.. 
The front free surface and rear free surface are in continual motion and can be tracked 
by the use of free-surface coordinates. The center-line boundary is fixed and its asso- 
ciated boundary conditions can be applied easily. 

Along the center-line boundary, boundary conditions can be applied by a fictitious 

column of nodal points on the opposite side of the longitudinal axis, each point being a 

mirror-image nodal point, and u^ t = 0. 

J, K 



Still boundary 


Figure 5.- Typical boundaries of target-projectile system during impact. 
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Boundary conditions along the free surfaces can be applied by assigning such prop- 
erty values to nodal points (adjacent to, but outside, the free surfaces) that the property 
values vary linearly from inside nodal points to outside nodal points and pass through the 
required free -surface values. 


Numerical Computation Cycle 

The dependent interior variables governing the fluid motions are u^j . , k , , , 

and p\ , and their governing equations are equations (7), (8), (9), (10), and (11). The 
J > K 

hydrostatic stress , is calculated from equation (13) and the equation of state. 

Jj k 

For the first approach, the dependent variables governing the motion of the free 
surfaces are 0^, u^ R , w^ R , and g, and their respective governing equations are 

equations (18), (23), and (24) and constant density along the free surfaces. 

9 p\ R 

The second approach is the same as the first approach except that -- is calcu- 
lated from equations (27), (28), (30), and (34). 


The stress derivatives are calculated from equations (25) and (26). 

When all the dependent variables i j}- . are specified at time t, the governing 

3 > K 

d Ak 

equations and boundary conditions are used to calculate — at time t. The values 
of the dependent variables at time t + r can then be found by use of the first two terms 
of a Taylor series: 



(35) 


Numerical Stability of the Free -Surface Coordinates 

A study of the stability of the equations governing the fluid interior motion can be 
found in the literature; for example, see reference 4. A complete study of the stability 
of the free-surface equations is beyond the scope of this report; however, when the free- 
surface equations are uncoupled from the interior flow (for the boundary treatment in the 
first approach) a study of their stability becomes tractable and is given in appendix C. 

NUMERICAL EXAMPLE 

Equations (7), (8), (9), (10), (11), (18), (23), and (24) were programed for the digital 
computer, and the results of a test case treated by the first approach are compared with 
the solution of a similar problem taken from reference 3. 
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The criteria used to evaluate the status and property values of nodal points exterior 
to the problem region are: 

(1) Nodal points closer to the free surfaces than Az/2 vertically or Ar/2 hori- 
zontally, as well as nodal points on the void side of the free surface, are considered to 
be exterior to the problem region. 

(2) Exterior nodal points adjacent to interior nodal points have such property values 
assigned to them that the boundary conditions are linearly satisfied. 

The test case simulates an aluminum truncated cone striking an aluminum target 
with an impact velocity of 20 x 10® cm/sec. The top diameter of the projectile is 
8 centimeters, the base diameter is 10 centimeters, and the height is 6 centimeters. The 
target is 2 centimeters thick. 

The case of reference 3 simulates an aluminum cylinder striking an aluminum tar- 
get with an impact velocity of 20 X 10® cm/sec. The cylinder diameter is 10 centimeters 
and the height is 10 centimeters. The target is 2 centimeters thick. 

The initial differences between the two problems are the sloping wall and shortened 
height of the test-case projectile. The sloping wall is necessary because infinite free- 
surface slopes cannot be accommodated by the equations derived herein. The difference 
in height of the two projectiles is inconsequential because at the time of comparison the 
rearward-traveling shock wave has not reached the rear free surface of the projectile. 

Figure 6 shows the projectile and target at impact. The direction of travel of the 
projectile is along the longitudinal coordinate into the target. 

Figure 7 shows the geometry of the two solutions 2.56 microseconds after impact. 
The dashed line represents the solution from reference 3 and the solid line represents 
the test-case solution. The most obvious difference between the solutions is the area 
inside the free surfaces. Since the solution from reference 3 accounts for the free sur- 
faces by diffusion of material into empty cells, it is to be expected that the free surfaces 
of the test case should be contained inside the free surfaces of reference 3, as is shown 
in figure 7. 

Figure 8 shows property variations along the center line. The velocity variations 
are shown in figure 8(a), density variations are shown in 8(b), and hydrostatic-stress 
variations are shown in 8(c). 

At the time of comparison, the shock wave traveling into the target has already 
encountered the target rear surface and the accompanying expansion wave has reentered 
the target and has overtaken the shock wave traveling into the projectile. The shock 
wave traveling into the projectile is situated at about z = 7 centimeters. 
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Reference 3 
Test case 



Figure 6.- Target-projectile system at impact. 








The velocity variations shown in figure 8(a) are in fair agreement. The mate- 
rial approaching the shock wave is being decelerated from its original speed of 
20 X 105 cm/sec. The material past the shock wave is being accelerated by the expan- 
sion wave from the rear free surface of the target. 

The density variations shown in figure 8(b) differ somewhat. The test-case density 
is lower (at the shock-wave location) than the density from reference 3 because of the 
value of H used in equation (13). Independent one -dimensional computations were per- 
formed to determine a value of H which would give smooth property variations , and a 
value of 5 was chosen. A smaller value of H would have resulted in a higher density 
peak for the test case. This dependence of property variations on H because of the 
small number of nodal points in the target has been indicated by several independent 
computer runs. For a smaller number of nodal points in the target, the dependence of 
property values on H is increased. The shape of the density curves near the free sur- 
face of the target differs appreciably. This difference is caused by the different methods 
of treating free surfaces . 

The stress curves shown in figure 8(c) are similar in shape but the solution of 
reference 3 is higher than the test-case solution by about lO 1 ^ dynes/cm2. This is 
thought to be caused by the value of H used in equation (13), as the peak stress, like 
the peak density, has been found to be a function of H. The stress curve for z less 
than 8 centimeters represents the shock wave traveling into the projectile and the stress 
curve for z greater than 8 centimeters represents the expansion wave overtaking the 
shock wave. The shock wave is about 2 centimeters thick in this case, far thicker than 
a real shock wave. This thickness is caused by the term q. A smaller value of H in 
equation (13) corresponds to thinner shock waves, and a larger value of H corresponds 
to thicker shock waves. 

The second approach was also programed (allowing for variable density along the 
free surfaces) and a computer run was made for the same initial conditions as were used 
with the first approach. The results with the second approach were nearly identical to 
those with the first approach except for density variations near the free surfaces. The 
different densities predicted by the two methods are presented in figure 9. 

Along the center line (fig. 9(a)) the density distributions are seen to be the same 
except near the rear free surface, where the second approach yields slightly lower values 
than the first approach. Figures 9(b) and 9(c) show density variations along the front 
and rear free surfaces, respectively. The curves from the second approach show small 
density increases at the positions where the radially traveling shock wave is interacting 
with the free surfaces, and small density decreases behind the shock wave. 

In choosing between the first and second approach, it should be noted that the shock- 
wave smearing and the adiabatic energy equation (no heat conduction) used herein are 
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responsible for the constant free -surface density in the first approach and for the small 
density variations along the free surfaces in the second approach. Since the energy equa- 
tion is adiabatic, no mechanism is available to bring heat to the free surface to cause 
expansion; furthermore, the shock-wave smearing decreases property gradients so 
drastically that no appreciable heat could be transferred to the free surfaces even if heat 
conduction terms were included in the energy equation. 

The first approach is mathematically correct for the working media considered 
and it is simpler to use than the second approach. 

The roughness in the density variation with the second approach for r between 4 
and 5 centimeters, shown in figure 9(b), is thought to be caused by the integration scheme 
used here. More computer runs have shown that this roughness can be eliminated by a 
small amount of smoothing of the data (a suitable weighting of the density at each point 
with its two nearest neighbors). 


CONCLUDING REMARKS 

A method of tracking free surfaces has been developed. Also, finite-difference 
approximations to the continuity equation and general transport equation have been derived 
in which the convective terms are simplified and mass, momentum, and internal energy 
are conserved in transit. The resulting fluid motions calculated by finite differences 
did not exhibit obvious instabilities. Two different methods are given for computing 
density variations along the free surfaces. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., August 27, 1969. 
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APPENDIX A 


FINITE-DIFFERENCE APPROXIMATION FOR THE TRANSPORT EQUATIONS 

IN CYLINDRICAL COORDINATES 

Consider the elemental volume 5j in figure 10, where 

n2ir 

6j = J rj Ar Az dO = 27rrj Ar Az 
and consider the density p\ , of a fluid flowing through 6i at 

Jj k J 

t = It 


r = jAr 
z = kAz 

The conservation of mass flowing through 8j requires that 

dp 1 . 


6j _iS + r wl/2) ^ i+(1/2);k ^ +(1/2)jk Az - 2* rj. (1/2) u'. (1/2) k p'j. (1/2)ik Az 


+ 2l r i w U + (l/2) Pi,k + (l/2) 4r ' 2,r r i w J,k-(l/2) p j,k-(l/2) ir * s (P)j,k 5 ) < Al > 



Radial coordinate, r 


Figure 10.- Elemental volume in cylindrical coordinates. 
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APPENDIX A 


where u is the radial velocity component, w is the longitudinal velocity component, 
and S(p) is a mass source per unit volume. From both sides of equation (Al) sub- 
tract the quantity 


J .1 


l l 


211 r j+(l/2) p j,k u j+(l/2),k Az - 2ir r j -(1/2) p j,k u j -(1/2) ,k Az 


+ 2 * r j 4k *i,k+(l/2) ar - 21 r i 4k 4k-(l/2) ar 


and after some rearrangement get 


6j _ik + 2, r. +(1/2) u I J+(1/2)>k (p? +(1/2)jk - P^a 


+ 2 » r j-(l/2) 4(1/2), k( P U - 4(1/2), k) 4z + 2 * r i 4 mi/2)(4mi/ 2) - p^ir 

+ 2w r s »i^-< i/2) ( p U - 4Mi/2)) ar = s( 4,k 6 i - 2 ’ r i 44^1/2),!, - 4(1/2), k) az 


U. 

- 21 r i 4kH,Ml/2) - ”i,Ml/2)) ar - 2 " p j,k 


Ar Az 


It is shown in appendix B that space differences, such as those on the left-hand side 
of equation (A2), should be taken in the opposite direction from the local velocity vector. 
After noticing that the conservation properties of equation (A2) are not changed by 
cycling all the first or second subscripts of jany term by the same amount, equation (A2) 
can be changed to 

dp 1 . 

6j + a 2 ^ rj u^ a/2) k (p' +a k - Pj )k ) az l t2 ”i 4Mb/2)(4k+b - 4/ r 
+ ” 4(a/2),k(4a,k - 4k) ar Az = S(p) i,k 6 i ' 2 ’ r i Pj,k(4(l/2),k ' 4(1/2), k) az 


2 ” r i 4k Ml/ 2 ) - w j,k-(l/2)) ar - 2 ' 4k 


i 4( 
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where a is sgn(-u) and b is sgn(-w), Equation (A3) was obtained from equation (A2) 
by cycling the j (or k) subscripts of the convective difference terms by ±1 so as to 
insure that each difference is taken in the opposite direction from the u (or w) velocity 
component. 

The effects of the index cycling can be seen by noticing (for example) that cycling 
the k-index on convective differences taken in the z-direction in equation (A2) is equiva- 
lent to adding the term 

-» r ) Ar [| w j ,k+(l/2) | (**j ,k+l - Pj,k)- | w i,k-(l/2)|(' J j,k - 

to equation (A2). When this term is associated with the time derivative of density in 
equation (A2) , there follows 



Thus the term added for the z-direction (to obtain eq. (A3) from eq. (A2)) acts like 
the second-order differential operator in the classical heat -conduction equation and helps 
smooth irregularities in the density variations in the z-direction. Finally, divide equa- 
tion (A3) by 5j and get 



By virtue of its derivation, equation (A4) conserves mass in transit. 

Consider next some quantity , which must be conserved in a time -integration 

1 J K 

scheme so that 
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which is identically equal to 


I( p j*k Q j|k) 6 i ‘ Q j,k) 6 j + Z( P Lk Q j,k) 5 i ' J( p j,k Q j,k) 5 j 

j,k j,k j,k j,k 


= Constant (A6) 


If the time derivative of , is taken to be 

P l> k _ p 3,k p ],k 

8t r 

then equations (A6) and (A7) can be combined to get 


8 ( P i,k Q j ,k ) . n i JSLk _o* 


at 


P 3 ?k — at” " V ~it~ " T ' 


at at 


= o 


(A7) 


(A8) 


which shows that the finite -difference approximation for the differentiation of products is 


p ],k at ^j,k at at at 


(A9) 


Qj k , which flows through 6j, to find with the aid of equations (A4) and (A9) that 


Finally, an approach similar to that used to derive equation (A4) can be applied to 


Q l . , - Q l . 


4 s ♦ • + b -Wi "W> 


vj, - 


Jy 


= s(pQ) i,k- r ^-^ (A10 > 


Equation (A10) conserves q { , in transit. 

J> K 
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RELATIONSHIP BETWEEN FINITE DIFFERENCES TAKEN AT 
CONSTANT TIME, AT CONSTANT SPACE POSITION, AND 
ALONG THE PARTICLE PATH 


Equation (4) shows the relationship between derivatives taken at constant time, 
derivatives taken at constant space position, and derivatives taken along the particle path. 
In order to establish a finite -difference approximation for equation (4), consider a particle 
having a describing property Q and moving along a particle trajectory. A time differ- 
ence taken at f is identically equal to 


Q(r , t+r) - Q(?,t) = [Q(r, t+r) - Q(?-/3, t)] + [Q(?-/3, t) - Q(r,tj] (Bl) 

so that 


Q(r, t+r) - Q(r-/3, t) = [Q(r, t+r) - Q(r,t)] + [Q(r,t) - Q(r-/3, t)] (B2) 


where 0 is an arbitrary vector. Equation (B2) is a finite -difference relationship equi- 
valent to equation (4). Thus the left-hand side of equation (B2) is a comoving difference. 
Since a particle trajectory which ends at (r, t+r) must begin at (r-Vr, t), where V is 
the local particle velocity, /3 must be equal to Vr. Denote the first difference in equa- 
tion (B2) by AQ tra . ectory , the second difference by AQ^ mg , and the third difference by 

AQgpace* Take the limit of equation (B2) as r goes to zero, and notice that 

AQ ^ 

usually is discontinuous at r. The value of — sgace mus t ^g calculated from the space 

m 

octant which contains the starting point of the particle trajectory that ends at (r, t+T). 

With this convention, equation (B2) becomes 


[AQ(r,t)] 

1 trajectory _ 1 

[AQ(r,t)] 

L aq| 

(time . 

K* 

) space 1 


T 

T 

+ 


/3 



(B3) 
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NUMERICAL STABILITY OF THE FREE -SURFACE EQUATIONS 

The numerical stability study of the free-surface equations is made tractable by 
holding o\ , constant in equations (25) and (26). Suppose u, w, and <fi are known 

numerical solutions for the free-surface motions. Let U^ , , and 4?^ be small space 

disturbances added to u^ B , Wj B , and in equations (18), (23), and (24). Neglect 
products of the small disturbance variables to find 


"S 41 - rf j ■ - a - «i) ■ - a - faW 


_L-_ L.B 

fa 8r 


' d d d 

*i , Vi ~ *1-1 

fa- fa 


^ - a SF U j)B ("W - ^ - » £(fa,B - 


d(j\ r, & 


r .i,B 

rl 3Z 


Z k-^J 


*i +1 »U«U - *j) - * &{fa - + "S 


where a is sgn(-u). Next let 


U Z j =^A Z n e inyr 

n 

Wj = ^ B^ e in rr 




e myr 
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l l 

iS. = a U j + - a ^ ."- U -l' g . 
9r Ar 


aw , »U.B - <B 

ar Ar 


<t> 1 . -<p l 

r t _I_Q I -t 


u = u, T, > 0 


(CIO) 


w = w. 


(Cll) 


1 9 ct 1 


P 8r n Z ar 

p j,B 


(C12) 


Substitute equations (C4) to (C12) into equations (Cl) to (C3) to get for each value 


of n 


A* +1 = A 1 - — u(l - e-^ArV - t fi V + 1 i2A_L_ + 2i , - sln ( n rfr) \ c i n (cl3) 

” n Ar v n n p8r U k -^ ^ + i-^.J n 


B i+1 = B* - 

n n Ar 


i-full. e- ln vAr) J , t |» a ! + t» _» c*„ 

n Ar v /D n ar n o az _ ^Z n 


ar n p a z _ , z 

Z k 


pZ+1 pZ T 

^n ~ 4 “ £7 




These equations can be written in matrix form as 


where 


#7 = [Mfr.TjK.Bk 
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m 12 = 0 




2i sin(nyAr)^ 


P9r k-^ 


M 2 i = -r f 


M. 


22 


= 1 ~ £u(l - e-^rAr) 


m 23 = J I 2 - 

4,5 P 9Z „ . 




“31 '- T ^ 


M 32 = t 

M 33 = 1 -A'‘(l- e ’ toyir ) 


and [M(n,T)J is called the amplification matrix. By repeated application of 
equation (C16), 


>1 N 

A n 


l+l 

ii 

Z+1 


r 

A i 


U[M(n,r)]W 


v. n ^ 






(C17) 


Now let r tend to zero in such a way that 

t > 0 
t = It 
1 - 0 0 

The integration scheme given by equation (C17) is said to be stable if and only if 


[M(n 


< Constant 


(CIS) 


for any positive constant and l sufficiently large, where Q±j is an arbitrary unit vec- 
tor and | | indicates the norm. 
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If jM(n, rjj can be written as 


[M(n,r)j g [M(n,0)j + [Constant] t 
as r tends to zero, and if the largest eigenvalue of jM(n,0)j 


(C19) 

is determined so that 


|x(m)|gl (m = 1, 2, 3) (C20) 

then (see ref. 4) equation (C18) holds true. Equation (C19) is seen to hold by inspection 
of equation (C16). Furthermore, the eigenvalues of [M(n,0)j are identical and equal to 

X(m) = l - £u(l -e- in rAr) 


Thus if 

J-ugi (C21) 

then fx(m)| g 1. Therefore equation (C17) is stable under numerical integration. 
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